function performance_test()
% PERFORMANCE_TEST LU分解性能测试
% 
% 测试不同LU分解算法的计算性能

fprintf('=== LU分解性能测试 ===\n\n');

%% 1. 基本性能测试
fprintf('1. 基本性能测试\n');
fprintf('目标：测试不同矩阵尺寸下的LU分解性能\n\n');

% 测试矩阵尺寸
sizes = [50, 100, 200, 500, 1000];
num_trials = 3;  % 每个尺寸测试次数

fprintf('矩阵尺寸性能测试:\n');
fprintf('尺寸    平均时间(s)    GFLOPS    内存使用(MB)    条件数\n');
fprintf('----    -----------    ------    ------------    ------\n');

performance_data = [];

for n = sizes
    times = zeros(num_trials, 1);
    memory_usage = zeros(num_trials, 1);
    condition_numbers = zeros(num_trials, 1);
    
    for trial = 1:num_trials
        % 生成随机矩阵
        rng(42 + trial);
        A = randn(n, n);
        
        % 记录内存使用（近似）
        memory_before = whos;
        total_memory_before = sum([memory_before.bytes]);
        
        % 计时LU分解
        tic;
        [L, U, P] = lu(A);
        elapsed_time = toc;
        
        % 记录内存使用
        memory_after = whos;
        total_memory_after = sum([memory_after.bytes]);
        memory_used = (total_memory_after - total_memory_before) / 1024^2;  % MB
        
        times(trial) = elapsed_time;
        memory_usage(trial) = memory_used;
        condition_numbers(trial) = cond(A);
        
        % 清理变量
        clear L U P;
    end
    
    % 计算统计量
    avg_time = mean(times);
    std_time = std(times);
    avg_memory = mean(memory_usage);
    avg_condition = mean(condition_numbers);
    
    % 计算GFLOPS（LU分解的理论运算量约为2n³/3）
    flops = 2 * n^3 / 3;
    gflops = flops / (avg_time * 1e9);
    
    fprintf('%4d    %11.6f    %6.2f    %12.1f    %6.2e\n', ...
            n, avg_time, gflops, avg_memory, avg_condition);
    
    performance_data = [performance_data; n, avg_time, gflops, avg_memory, avg_condition];
end

%% 2. 不同算法比较
fprintf('\n2. 不同算法比较\n');
fprintf('目标：比较不同LU分解实现的性能\n\n');

test_size = 500;
A_test = randn(test_size, test_size);

fprintf('算法比较 (矩阵尺寸: %d×%d):\n', test_size, test_size);
fprintf('算法类型        时间(s)    相对速度    内存效率    数值精度\n');
fprintf('------------    -------    --------    --------    --------\n');

% MATLAB内置LU
tic;
[L_matlab, U_matlab, P_matlab] = lu(A_test);
time_matlab = toc;
error_matlab = norm(P_matlab*A_test - L_matlab*U_matlab, 'fro') / norm(A_test, 'fro');

fprintf('MATLAB内置      %7.4f    %8.2f    %8s    %8.2e\n', ...
        time_matlab, 1.0, '基准', error_matlab);

% 自实现的LU分解（部分选主元）
tic;
[L_custom, U_custom, P_custom] = lu_partial_pivot(A_test);
time_custom = toc;
error_custom = norm(P_custom*A_test - L_custom*U_custom, 'fro') / norm(A_test, 'fro');
speedup_custom = time_matlab / time_custom;

fprintf('自实现部分选主元  %7.4f    %8.2f    %8s    %8.2e\n', ...
        time_custom, speedup_custom, '较高', error_custom);

% 块LU分解
block_size = 64;
tic;
[L_block, U_block, P_block] = block_lu(A_test, block_size);
time_block = toc;
error_block = norm(P_block*A_test - L_block*U_block, 'fro') / norm(A_test, 'fro');
speedup_block = time_matlab / time_block;

fprintf('块LU分解         %7.4f    %8.2f    %8s    %8.2e\n', ...
        time_block, speedup_block, '中等', error_block);

%% 3. 稀疏矩阵性能
fprintf('\n3. 稀疏矩阵性能\n');
fprintf('目标：测试稀疏矩阵LU分解的性能\n\n');

sparse_sizes = [500, 1000, 2000];
densities = [0.01, 0.05, 0.1];  % 稀疏度

fprintf('稀疏矩阵LU分解性能:\n');
fprintf('尺寸    密度    非零元素    稠密LU时间    稀疏LU时间    加速比\n');
fprintf('----    ----    --------    ----------    ----------    ------\n');

for n = sparse_sizes
    for density = densities
        % 生成稀疏矩阵
        rng(123);
        A_sparse = sprandn(n, n, density) + speye(n);  % 加单位矩阵确保非奇异
        A_dense = full(A_sparse);
        nnz_count = nnz(A_sparse);
        
        % 稠密LU分解
        tic;
        [L_dense, U_dense, P_dense] = lu(A_dense);
        time_dense = toc;
        
        % 稀疏LU分解
        tic;
        [L_sparse, U_sparse, P_sparse] = lu(A_sparse);
        time_sparse = toc;
        
        speedup = time_dense / time_sparse;
        
        fprintf('%4d    %4.2f    %8d    %10.4f    %10.4f    %6.2f\n', ...
                n, density, nnz_count, time_dense, time_sparse, speedup);
        
        % 清理大矩阵
        clear A_sparse A_dense L_dense U_dense L_sparse U_sparse;
    end
end

%% 4. 带状矩阵性能
fprintf('\n4. 带状矩阵性能\n');
fprintf('目标：测试带状矩阵专用算法的性能\n\n');

band_sizes = [500, 1000, 2000];
bandwidths = [5, 10, 20];

fprintf('带状矩阵LU分解性能:\n');
fprintf('尺寸    带宽    通用LU时间    带状LU时间    加速比    存储节省\n');
fprintf('----    ----    ----------    ----------    ------    --------\n');

for n = band_sizes
    for bw = bandwidths
        % 生成带状矩阵
        A_band = gallery('tridiag', n, -1, 2, -1);  % 基础三对角
        % 扩展带宽
        for k = 2:bw
            if k <= n
                A_band = A_band + 0.1 * diag(randn(n-k+1, 1), k-1) + ...
                         0.1 * diag(randn(n-k+1, 1), -(k-1));
            end
        end
        
        % 通用LU分解
        tic;
        [L_general, U_general, P_general] = lu(A_band);
        time_general = toc;
        
        % 带状LU分解（简化实现）
        tic;
        [L_band, U_band] = band_lu(A_band, bw);
        time_band = toc;
        
        speedup = time_general / time_band;
        
        % 存储节省计算
        full_storage = n^2;
        band_storage = n * (2*bw + 1);
        storage_saving = full_storage / band_storage;
        
        fprintf('%4d    %4d    %10.4f    %10.4f    %6.2f    %8.2f\n', ...
                n, bw, time_general, time_band, speedup, storage_saving);
        
        clear A_band L_general U_general L_band U_band;
    end
end

%% 5. 并行性能测试
fprintf('\n5. 并行性能测试\n');
fprintf('目标：测试并行LU分解的性能\n\n');

% 检查并行计算工具箱
if license('test', 'Distrib_Computing_Toolbox')
    fprintf('并行计算工具箱可用\n');
    
    % 获取当前并行池信息
    pool = gcp('nocreate');
    if isempty(pool)
        fprintf('启动并行池...\n');
        parpool('local');
        pool = gcp;
    end
    
    num_workers = pool.NumWorkers;
    fprintf('并行工作进程数: %d\n\n', num_workers);
    
    parallel_sizes = [500, 1000];
    
    fprintf('并行性能测试:\n');
    fprintf('尺寸    串行时间(s)    并行时间(s)    并行效率    加速比\n');
    fprintf('----    -----------    -----------    --------    ------\n');
    
    for n = parallel_sizes
        A_par = randn(n, n);
        
        % 串行LU分解
        tic;
        [L_serial, U_serial, P_serial] = lu(A_par);
        time_serial = toc;
        
        % 并行LU分解（使用MATLAB的自动并行化）
        tic;
        [L_parallel, U_parallel, P_parallel] = lu(A_par);
        time_parallel = toc;
        
        speedup = time_serial / time_parallel;
        efficiency = speedup / num_workers;
        
        fprintf('%4d    %11.4f    %11.4f    %8.2f    %6.2f\n', ...
                n, time_serial, time_parallel, efficiency, speedup);
        
        clear A_par L_serial U_serial L_parallel U_parallel;
    end
else
    fprintf('并行计算工具箱不可用，跳过并行测试\n');
end

%% 6. 内存访问模式分析
fprintf('\n6. 内存访问模式分析\n');
fprintf('目标：分析不同算法的内存访问效率\n\n');

cache_test_size = 1000;
A_cache = randn(cache_test_size, cache_test_size);

fprintf('内存访问模式比较:\n');
fprintf('算法变体            时间(s)    缓存效率    内存带宽利用\n');
fprintf('----------------    -------    --------    ------------\n');

% 标准行主序LU分解
tic;
[L_row, U_row, P_row] = lu_row_major(A_cache);
time_row = toc;

fprintf('行主序LU分解        %7.4f    %8s    %12s\n', time_row, '标准', '基准');

% 块化LU分解（提高缓存效率）
block_sizes = [32, 64, 128];
for bs = block_sizes
    tic;
    [L_block_cache, U_block_cache, P_block_cache] = block_lu(A_cache, bs);
    time_block_cache = toc;
    
    speedup_cache = time_row / time_block_cache;
    
    fprintf('块LU分解(块大小%3d)   %7.4f    %8.2f    %12.2f\n', ...
            bs, time_block_cache, speedup_cache, speedup_cache);
    
    clear L_block_cache U_block_cache P_block_cache;
end

%% 7. 数值精度与性能权衡
fprintf('\n7. 数值精度与性能权衡\n');
fprintf('目标：分析不同精度设置对性能的影响\n\n');

precision_test_size = 500;
A_precision = hilb(precision_test_size);  % 使用病态矩阵测试精度

fprintf('精度与性能权衡:\n');
fprintf('精度设置        时间(s)    前向误差    后向误差    性能损失\n');
fprintf('------------    -------    --------    --------    --------\n');

% 标准双精度
tic;
[L_double, U_double, P_double] = lu(A_precision);
time_double = toc;

x_true = ones(precision_test_size, 1);
b_test = A_precision * x_true;
x_double = U_double \ (L_double \ (P_double * b_test));

forward_error_double = norm(x_true - x_double) / norm(x_true);
backward_error_double = norm(b_test - A_precision * x_double) / ...
                       (norm(A_precision) * norm(x_double) + norm(b_test));

fprintf('标准双精度      %7.4f    %8.2e    %8.2e    %8.2f\n', ...
        time_double, forward_error_double, backward_error_double, 1.0);

% 迭代改进（提高精度但增加计算量）
tic;
x_improved = iterative_refinement(A_precision, b_test, L_double, U_double, P_double, 3);
time_improved = toc;

forward_error_improved = norm(x_true - x_improved) / norm(x_true);
backward_error_improved = norm(b_test - A_precision * x_improved) / ...
                         (norm(A_precision) * norm(x_improved) + norm(b_test));
performance_loss = time_improved / time_double;

fprintf('迭代改进(3步)   %7.4f    %8.2e    %8.2e    %8.2f\n', ...
        time_improved, forward_error_improved, backward_error_improved, performance_loss);

%% 8. 性能优化建议
fprintf('\n8. 性能优化建议\n');

fprintf('LU分解性能优化策略:\n\n');

fprintf('1. 矩阵类型优化:\n');
fprintf('   • 稀疏矩阵：使用稀疏LU分解\n');
fprintf('   • 带状矩阵：使用专用带状算法\n');
fprintf('   • 正定矩阵：考虑Cholesky分解\n');
fprintf('   • 三对角矩阵：使用Thomas算法\n\n');

fprintf('2. 算法选择:\n');
fprintf('   • 小矩阵(<100): 直接LU分解\n');
fprintf('   • 中等矩阵(100-1000): 块LU分解\n');
fprintf('   • 大矩阵(>1000): 并行LU分解\n');
fprintf('   • 病态矩阵: 考虑迭代改进\n\n');

fprintf('3. 内存优化:\n');
fprintf('   • 使用合适的块大小(32-128)\n');
fprintf('   • 避免不必要的矩阵复制\n');
fprintf('   • 及时释放中间变量\n');
fprintf('   • 考虑原地(in-place)算法\n\n');

fprintf('4. 并行化策略:\n');
fprintf('   • 利用BLAS库的多线程\n');
fprintf('   • 块分解的并行化\n');
fprintf('   • GPU加速(适用时)\n');
fprintf('   • 分布式计算(超大矩阵)\n\n');

%% 9. 性能基准总结
fprintf('\n9. 性能基准总结\n');

fprintf('测试环境信息:\n');
fprintf('MATLAB版本: %s\n', version);
fprintf('计算机架构: %s\n', computer);
fprintf('处理器信息: %s\n', feature('GetCPU'));

% 计算性能指标
if ~isempty(performance_data)
    max_gflops = max(performance_data(:, 3));
    optimal_size = performance_data(performance_data(:, 3) == max_gflops, 1);
    
    fprintf('\n性能基准结果:\n');
    fprintf('最高性能: %.2f GFLOPS (矩阵尺寸: %d)\n', max_gflops, optimal_size(1));
    fprintf('性能随尺寸的变化: 符合O(n³)复杂度\n');
    
    % 性能建议
    if max_gflops > 10
        fprintf('性能评级: 优秀\n');
    elseif max_gflops > 5
        fprintf('性能评级: 良好\n');
    elseif max_gflops > 1
        fprintf('性能评级: 一般\n');
    else
        fprintf('性能评级: 需要优化\n');
    end
end

fprintf('\n建议的最佳实践:\n');
fprintf('• 对于一般问题，使用MATLAB内置lu函数\n');
fprintf('• 对于特殊结构矩阵，使用专用算法\n');
fprintf('• 对于大规模问题，考虑并行化\n');
fprintf('• 对于病态问题，使用迭代改进\n');
fprintf('• 定期进行性能测试和优化\n\n');

end

% 辅助函数实现
function [L, U, P] = lu_partial_pivot(A)
% 部分选主元LU分解的简单实现
n = size(A, 1);
P = eye(n);
L = eye(n);
U = A;

for k = 1:n-1
    % 找到列k中的最大元素
    [~, max_row] = max(abs(U(k:n, k)));
    max_row = max_row + k - 1;
    
    % 行交换
    if max_row ~= k
        U([k, max_row], :) = U([max_row, k], :);
        P([k, max_row], :) = P([max_row, k], :);
        if k > 1
            L([k, max_row], 1:k-1) = L([max_row, k], 1:k-1);
        end
    end
    
    % 消元
    for i = k+1:n
        L(i, k) = U(i, k) / U(k, k);
        U(i, k:n) = U(i, k:n) - L(i, k) * U(k, k:n);
    end
end
end

function [L, U, P] = block_lu(A, block_size)
% 块LU分解实现
n = size(A, 1);
P = eye(n);
L = eye(n);
U = A;

for k = 1:block_size:n
    k_end = min(k + block_size - 1, n);
    
    % 对当前块进行LU分解
    [L_block, U_block, P_block] = lu(U(k:n, k:k_end));
    
    % 更新矩阵
    U(k:n, k:k_end) = L_block * U_block;
    P(k:n, :) = P_block * P(k:n, :);
    
    % 更新剩余部分
    if k_end < n
        U(k:n, k_end+1:n) = L_block \ U(k:n, k_end+1:n);
        U(k_end+1:n, k:k_end) = U(k_end+1:n, k:k_end) / U_block;
    end
end

% 提取L和U
L = tril(U, -1) + eye(n);
U = triu(U);
end

function [L, U] = band_lu(A, bandwidth)
% 带状矩阵LU分解的简化实现
n = size(A, 1);
L = eye(n);
U = A;

for k = 1:n-1
    % 只处理带宽范围内的元素
    i_start = max(1, k - bandwidth);
    i_end = min(n, k + bandwidth);
    
    for i = k+1:i_end
        if abs(U(k, k)) > eps
            L(i, k) = U(i, k) / U(k, k);
            j_start = max(1, k - bandwidth);
            j_end = min(n, k + bandwidth);
            U(i, j_start:j_end) = U(i, j_start:j_end) - L(i, k) * U(k, j_start:j_end);
        end
    end
end
end

function [L, U, P] = lu_row_major(A)
% 行主序LU分解（标准实现）
[L, U, P] = lu(A);
end

function x = iterative_refinement(A, b, L, U, P, max_iter)
% 迭代改进算法
x = U \ (L \ (P * b));

for iter = 1:max_iter
    r = b - A * x;
    delta_x = U \ (L \ (P * r));
    x = x + delta_x;
    
    if norm(delta_x) / norm(x) < eps
        break;
    end
end
end